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ABSTRACT 

We present a method for calculating the infrared emission from a population of dust 
grains heated by starlight, including very small grains for which stochastic heating by 
starlight photons results in high temperature transients. Because state-to-state tran- 
sition rates are generally unavailable for complex molecules, we consider model PAH, 
graphitic, and silicate grains with realistic vibrational mode spectra and realistic ra- 
diative properties. The vibrational density of states is used in a statistical-mechanical 
description of the emission process. Unlike previous treatments, our approach fully in- 
corporates multiphoton heating effects, important for large grains or strong radiation 
fields. 

We discuss how the "temperature" of the grain is related to its vibrational energy. By 
comparing with an "exact" statistical calculation of the emission process, we determine 
the conditions under which the "thermal" and the "continuous cooling" approximations 
can be used to calculate the emission spectrum. 

We present results for the infrared emission spectra of PAH grains of various sizes 
heated by starlight. We show how the relative strengths of the 6.2, 7.7, and 11.3/im 
features depend on grain size, starlight spectrum and intensity, and grain charging 
conditions. We show results for grains in the "cold neutral medium", "warm ionized 
medium", and representative conditions in photodissociation regions. Our model results 
are compared to observed ratios of emission features for the Milky Way and other 
galaxies, and for the M17 and NGC 7023 photodissociation regions. 

Subject headings: Galaxies: ISM; Infrared: ISM: Continuum; ISM: Dust 

1. Introduction 



Very small grains - consisting of tens to hundreds of atoms - play a major role in the infrared 
emission from the interstellar medium. These grains are small enough that the time-averaged vibra- 
tional energy (E) is smaller than or comparable to the energy of the starlight photons which heat 
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the grains. Stochastic heating by absorption of starhght therefore results in transient "temperature 
spikes", during which much of the energy deposited by the starhght photon is reradiated in the 
infrared. 

The idea of transient heating of very small grains was first introduced by Greenberg (1968). 
Its importance was not apparent until the detection of the near infrared emission of reflection 
nebulae (Sellgren et al. 1983) and detection by the Infrared Astronomical Satellite (IRAS) of 
12 and 25/im Galactic emission which was far in excess of the emission expected for interstellar 
dust at T 20K (Boulanger & Perault 1988). Subsequent measurements by the Diffuse Infrared 
Background Experiment (DIRBE) instrument on the Cosmic Background Explorer (COBE) satellite 
confirmed this and detected additional broadband emission at 3.5 and 4.9/xm (Arendt et al. 1998). 
More recently, spectrometers aboard the Infrared Telescope in Space (IRTS) (Onaka et al. 1996; 
Tanaka et al. 1996) and the Infrared Space Observatory (ISO) (Mattila et al. 1996) have shown 
that the diffuse interstellar medium radiates strongly in emission features at 3.3, 6.2, 7.7, 8.6, and 
11.3)um. These emission features, previously observed from a small number of bright reflection 
nebulae, planetary nebulae, and HII regions, have been tentatively identified as emission from 
polycyclic aromatic hydrocarbons, or PAHs (Leger k, Puget 1984; AUamandola, Tielens, k, Barker 
1985). 

To understand these observations we need to be able to calculate the infrared emission expected 
for a given model of interstellar dust. The question we address here is: what approximations can be 
used to carry out accurate calculations of thermal emission from a population of very small grains 
exposed to stochastic heating by starlight? 

There have been a number of previous studies of the discrete heating of very small grains (Duley 
1973; Greenberg & Hong 1974; Gail & Seldmayr 1975; Greenberg 1976; Purcell 1976; Drapatz & 
Michel 1977; Andriesse 1978; Aannestad & Kenyon 1979; Leger & Puget 1984; Draine & Anderson 
1985; Puget et al. 1985; Desert et al. 1986; Dwek 1986; Guhathakurta & Draine 1989; Aannestad 
1989; Siebenmorgen et al. 1992; Manske & Henning 1998; Duley & Poole 1998). Most studies 
have made three approximations: (1) the "thermal" approximation for excitation; (2) use of "bulk" 
specific heats; and (3) the "continuous cooling" approximation. While accurate for large grains, 
these approximations could potentially introduce significant inaccuracies for very small grains. 

Barker & Cherchneff (1989) and d'Hendecourt et al. (1989) found that the thermal approxi- 
mation was valid for computing thermal emission spectra. AUamandola, Tielens, & Barker (1989), 
Schutte, Tielens, & AUamandola (1993), and Cook & Saykally (1998) compared calculations of 
infrared emission from PAHs using the thermal approximation with the results of a more detailed 
statistical treatment; they found that the thermal approximation provided a good approximation 
for the emission from low- frequency modes, although it overestimated the emission from high fre- 
quency modes such as the 3.3/im C-H stretching mode. All of these studies used the "continuous 
cooling" approximation to discuss the overall emission from a PAH following single-photon heating. 

The validity of the "continuous cooling" approximation has been questioned by Siebenmorgen 
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et al. (1992), Manske & Henning (1998), and Duley & Poole (1998). Siebenmorgen et al. (1992) 
and Manske & Henning (1998) attempted to treat the cooling process as discrete events, but did 
not correctly evaluate the transition probabilities, as we show below in §6.1. 

The present paper addresses the vibrational excitation of interstellar grains, both small and 
large, and the resulting infrared emission. If energy levels and transition probabilities were known, 
we could (at least in principle) solve for the statistical steady-state populations of the different 
energy levels of grains illuminated by a known radiation field. However, this level of detailed infor- 
mation is generally unavailable, for even the smallest and simplest PAH molecules. We therefore 
construct a "model PAH" grain, C„H^, with a realistic spectrum of energy levels. Following pre- 
vious workers (e.g., Schutte et al. 1993, Boulanger et al. 1998) we adopt realistic absorption cross 
section Cabs (A), depending on the number of C atoms, the H/C ratio, and whether the PAH is neu- 
tral or ionized. We construct a similar model for silicate grains, with appropriate mode spectrum 
and absorption cross section. We divide the energy level spectrum into a manageable number of 
"bins", estimate the "bin-to-bin" transition probabilities, and can then solve for the steady-state 
probability Pj of occupying energy bin j. 

Our analysis allows for discrete quantum states at low vibrational energies, and what is ef- 
fectively a continuum of vibrational levels at large energies. As we show, the only simplifying 

assumptions that need to be made are: (1) the absorption cross section Cabs (A) is independent 
of the degree of vibrational excitation, and (2) internal vibrational relaxation rapidly redistributes 
internal energy among the vibrational degrees of freedom of the grain. While not exact, these are 
both expected to be good approximations. Once Cabs(-^) is adopted and the spectrum of vibra- 
tional "normal modes" (and from this the vibrational density of states) is specified, the physics 
of absorption and emission of radiation by the grain is essentially fully determined. With these 
assumptions, we can solve for the energy distribution function P{E) for a given grain in an arbi- 
trary radiation field, with full inclusion of "multiphoton heating" effects which become important 
for large grains or intense radiation fields. Because we solve for P{E) including all radiative tran- 
sitions, we obtain the "exact" solution against which we can then test solutions obtained using the 
thermal approximation and the continous cooling approximation. 

In §2 we estimate the spectrum of vibrational modes for small PAH and silicate grains, and 
obtain the vibrational density of states for such grains. The transition probabilities depend on the 
photon absorption cross sections, and in §3 we describe the cross sections adopted for PAH grains. 

We formulate the "exact" statistical treatment in §4. In §5 we discuss how the "temperature" 
of a grain is related to its vibrational energy content, and in §6 we show how the "thermal" 
approximation is related to the exact treatment. The radiative cooling time for vibrationally- 
excited grains is estimated as a function of E and grain size in §7. 

In §8 we discuss the method used to solve for the steady-state energy distribution. In §9 we 
show how the infrared emission spectrum can be calculated from the energy distribution, both 
for the exact statistical treatment and using the thermal approximation. We show some sample 
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solutions in §10, where we assess the accuracy of the "thermal" approximation. We find the 

thermal approximation to be quite accurate for calculation of the overall emission spectrum from a 
population of very small grains. Furthermore, we find that the "continuous cooling" approximation 
is itself quite accurate, and can therefore be used to greatly reduce the amount of computation 
required to find the temperature distribution functions. 

The above studies therefore provide a methodology for computing the infrared emission from 
dust grains irradiated by starlight. The method is accurate, and can be used for any specified 
mixture of grain sizes and compositions; we present results for both carbonaceous and silicate 
materials. We discuss our results in §11 and summarize the main points in §12. 

The paper is largely concerned with the physics of infrared emission from stochastically-heated 
grains, and practical methods for calculating the emission spectrum. Readers who are interested 
only in the resulting emission spectrum may wish to proceed directly to §10. In a following paper 
(Li & Drainc 2001a) wc apply this methodology to construct models to account for the observed 
infrared emission from the diffuse interstellar medium. 



2. Vibrational Mode Spectrum and Level Degeneracies 

A grain containing Na atoms will have 3 translational degrees of freedom, 3 rotational degrees 
of freedom, and iV^ = 3Na — 6 vibrational modes. Each vibrational mode j = I, •••,iVm has a 
fundamental frequency Uj. If we approximate the vibrational modes as harmonic oscillators, then 
the vibrational energy of the grain (not including zero-point energy) is 

Nm 

E = Y^ Vjhwj , (1) 

where the integer Vj > is the vibrational quantum number for mode j. If the vibrational system 
is in thermal equilibrium at temperature T, the expectation value of the energy E is 



where k is Boltzmann's constant. The heat capacity 



dE 



N„ 



e 



hwj/kT 



dT ^ [ 1 - exp (- /kT) 



hujj/kT 



(3) 



Ideally, we would have accurate knowledge of the N^n frequencies loj. Since the complete mode 
spectrum is known only for a small number of grain candidates such as C24H12 coronene (see 
below), we need to estimate the normal mode spectrum for candidate grains of different sizes and 
materials. 
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2.1. Normal Modes: Polycyclic Aromatic Hydrocarbons 

A polycyclic aromatic hydrocarbon molecule containing A^c C atoms and A'^h H atoms has 
= ^{Nh + Nq — 2) distinct vibrational modes. These can be separated into 3(iVc — 2) modes of 
the C-C skeleton, and SiVn modes associated with the C-H bonds at the periphery. The frequencies 
of these normal modes have been computed for a small number of polycyclic aromatic hydrocarbons, 
with some frequencies determined experimentally, but mode spectra are not yet available for most 
PAHs of interest. Here we estimate the mode spectrum for "generic" PAHs. We will treat 5 
different types of vibration separately: out-of-plane C-C modes, in-plane C-C modes, out-of-plane 
C-H bending, in-plane C-H bending, and C-H stretching. 

Graphite is a useful guide to the C-C modes. Krumhansl &: Brooks (1953) found that the lattice 
vibration specific heat of graphite could be approximated by a model where the out-of-plane and 
in-plane C-C vibrations each had a mode spectrum given by a two-dimensional Debye model, with 
Debye temperatures Qop ~ 950 K for the out-of-plane modes, and Qip ~ 2500 K for the in-plane 
modes. This suggests that the C-C vibrational modes for a PAH might be similarly approximated. 

A n-dimensional Debye spectrum has modes uniformly distributed in E'" up to a maximum 
energy A;G. Consider a mode spectrum of the form 



(1 - f3n] 



m 



1/n 



i = i,-,iv^. (4) 



For the PAH C-C modes we set the dimensionality n = 2, with = A'^c — 2 for the out-of-plane 
modes, and = 2(Arc — 2) for the in-plane modes. 

If Sj is independent of j, the modes are uniformly distributed in E"^. We will in fact take 

Sj = 1/2 for jV2or3 (5) 
= 1 for j = 2 or 3 (6) 

because this shift of modes 2 and 3 to lower frequency improves the match to the actual mode 
spectrum of coronene. 

As discussed in Appendix A we assume PAHs to be planar for Nq < 54, approximately spherical 
for A^c > 102, with intermediate shape for 54 < Nq < 102. Thus for A^c < 54 we expect the lowest 

— 1/2 . 

frequency mode wi oc , approxii 

A^c > 102. We achieve this by taking 



— \/2 1/3 

frequency mode wi oc A^,^ , approximately constant wi for 54 < A^c < 102, and ivi oc for 



^2 = for Nc < 54 

1 fNc - 54 



(2A^^ - 1) V 52 
1 



(2Ar^ - 1) 



(A^c - 2) fl02y^^ ^ 



52 V^cy 



for 54 < ATc < 102 

for A^c > 102 (7) 
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We use equations (4-7) with = Nq — 2 and kQop/hc = 600 cm~^ {Qop = 863 K) for the out-of- 
planc C-C modes, and = 2{Nc - 2) and k@ip/hc = 1740 cm~^ {@ip = 2504 K) for the in-plane 
C-C modes. 

We can adequately approximate the C-H stretching and bending modes by A'h out-of-plane 
bending modes at A^^ = (11.3;um)~^ = 886 cm~^, A'^h in-plane bending modes at A^^^^ = 
(8.6/im)~^ = 1161 cm~^, and A^h C-H stretching modes at ^clistr ~ (3.3/xm)~^ = 3030cm~^. 

For our model "astronomical PAH" sequence, we will assume 

A^H = int (0.57Vc + 0.5) Nc < 25 

= int (2.5^/iVc + 0.5) 25 < iVc < 100 

= int (0.25iVc + 0.5) A^c > 100 (8) 

where int(.x) is the integer part of x. The pericondensed species coronene C24H12; circumcoronene 
C54H18, and dicircumcoronene C96H24 are members of the sequence prescribed by eq. (8). Beyond 
A'c = 100 we assume H/C « 0.25. For a carbonaceous grain containing Nc carbon atoms, we 
associate a representative "radius" a = 10 A(A/c/468)^/^, corresponding to the carbon density 
p = 2.24 gcm~^ of graphite. 

We use C24H12 coronene, for which the mode spectrum is known (Cyvin 1982; Cyvin et al. 
1984), to test the above model. As shown in Figure 1, our "synthetic" mode spectrum is in excellent 
agreement with the actual mode spectrum of coronene. 

Applying our synthetic mode spectrum for PAHs in the limit Nc 00 and Nu = 0, the mode 
spectrum of the C-C modes reduces to the 2-dimensional Debye model. The heat capacity becomes 

Cgraph = (iVc - 2)k [/^(r/863K) + 2/^(r/2504K)] (9) 

U{x) ^- f ff\ , ax) = ^Ux) . (10) 

n Jo exp(y/a;) — 1 ax 

In Figure 2 we compare eq. (9) with experimental data for (bulk) graphite. The excellent agreement 
confirms that our model spectrum for the C-C vibrational modes applies from small PAH molecules 
up through macroscopic samples of graphite. 



2.2. Normal Modes: Silicate Grains 

In the case of silicate grains, we use experimental specific heats as a guide to the vibrational 
mode spectrum. The specifics heat of Si02 glass, obsidian glass (75% Si02 and 25% metal oxides 
by mass), and basalt glass (50% Si02, 50% metal oxides by mass) reported by Leger, Jura, & 
Omont (1985) are plotted Figure 2. As seen in Figure 2, the experimental results for bulk silicates 
(basalt and obsidian) can be satisfactorily reproduced if 2/3 of the vibrational modes are distributed 
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Fig. 1. — Njnodcs{E), the number of vibrational modes with energy < E, for coronene C24H12 (solid lines) compared 
to the model spectrum (broken lines) given by eq.(4-7) with A*'™ — 22 modes with k@op/hc = 600 cm^^, Nm = 44 
modes with kQip/hc — 1740 cm^ Nm = 12 modes with A"^ = 886, 1161, and 3030 cm (see text). 



according to a Debye model with n = 2 and Debye temperature = 500 K, and 1/3 of the modes 
are described by a Debye model with n = 3 and @ = 1500 K: 

Csii = {Na - 2)k [2/^(r/500K) + /^(r/1500K)] (11) 

where Na is the number of atoms in the cluster, and fn is given by eq. (10). Eq. (11) agrees quite 
well with the measured specific heats for obsidian and basalt, so we use this model to estimate the 
mode spectrum for silicate grains. 

— 1/3 

We expect the lowest frequency mode loi oc Na ■ To achieve this scaling for modes distributed 
according to an n— dimensional Debye model [eq. (4)], we set 

iV^n 1 

2iV^-l ' ^ ^ 

note that for n = 3 this gives P3 = 0. 
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Fig. 2. — Heat capacity of graphite (deSorbo & Tyler 1953; Furukawa & Douglas 1972) and model specific heat 
given by eq. (9); and heat capacity of Si02 glass, obsidian glass, and basalt glass (Leger, Jura, & Omont 1985), 
together with the model specific heat given by eq. (11). 



For small silicate grains, two-thirds of the vibrational modes are then approximated by eqs. (4- 
6) with n = 2, Qd = 500 K, and Nm = 2{Na — 2), and the remaining third by eqs. (4-6) with n = 3, 
@D = 1500 K, and Nm = Na — 2. The resulting normal mode spectrum is shown in Figure 3. 



2.3. Density of States 

Let there be a total of = 3{Na — 2) distinct normal modes. The general vibrational state can 
then be specified by the A^^-tuple (^1,^2, giving the vibrational quantum numbers for each 

of these modes. If each mode is approximated by a harmonic oscillator, then N{E), the number of 
distinct vibrational states with total vibrational energy less than or equal to E, can be calculated 
using the Beyer-Swinehart algorithm (Beyer & Swinehart 1973; Stein & Rabinovitch 1973). In 
Figure 4 we show N{E) computed for C24H12, using both the actual normal mode spectrum for 
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Fig. 3. — Vibrational mode spectrum adopted for silicate clusters with Na ~ 20, 35, and 100 atoms. 

coronene and our model normal mode spectrum for C24H12 (see Figure 1). The resulting N{E) are 
essentially identical for E/hc > 300 cm^^. 

The Beyer-Swinehart algorithm is remarkably efficient, and with IEEE standard 64-bit arith- 
metic can be applied to calculate N{E) up to E/hc = 10^ cm~^ for clusters with up to ~600 atoms. 
For a given cluster, we employ the Beyer-Swinehart algorithm for energies up to an energy Et 
where N{Et) « 10^'^'', beyond which we are limited by inability to calculate floating point numbers 
exceeding 2^^'^'^ ~ lO^'^^. For E > Et we calculate the density of states by integrating dS = dQ/T: 

lniV(£)=l„iV(£,)+y^_— , (13) 

where T{E) is the temperature at which the grain has energy E (see §5). Figures 4 and 5 show 
N{E) evaluated for C24H12, C400H100, C800H20O) and C4000H1000 using our model density of states. 
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3. Absorption Cross Sections 

Following previous workers (e.g., Schutte et al. 1993), we estimate the absorption cross sections 
for both neutral and ionized PAHs from far ultraviolet to far infrared (see Li & Draine 2001a for 
details) based on available experimental data (Allamandola et al. 1999; Hudgins &: Allamandola 
1999; and references therein) and guided by astronomical observations (e.g., Boulanger et al. 1998). 
In the ultraviolet and infrared the resulting cross sections are mainly a collection of Drude profiles: 
the a — a* transition mode (A~^ ~ 14//m~^); the tt — tt* transition (A~^ ~ 4.6/im~^); the C- 
H stretching mode (A = 3.3//m); the C-C stretching modes (A = 6.2,7.7/im); the C-H in-plane 
bending mode (A = 8.6/im); the C-H out-of-plane bending modes (A = 11.3, 11.9, 12. 7;um); and a 
few weak features attributed to the C-C bending modes (A = 16.4, 18.3, 21.2, 23.1/im). In addition 
to these features, we have included continuous absorption in the far-UV, near-UV/ visible, as well as 
a small amount of continuous absorption in the infrared. [For further information see Li &: Draine 
(2001)]. In Figure 6 we display the adopted absorption cross sections for neutral and ionized PAHs. 
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Fig. 5. — Total number of vibrational states N{E) with energy less than E for C400H100 and C8ooH2oo- Solid line 
is the direct Beyer-Swinehart calculation for A'^ < 10^""; broken line is from eq. (13). 



Following Li & Draine (2001), we make a smooth transition from PAH optical properties for 
^ 6 X 10^ (a < 50 A) to graphite properties for iVc > 5 x 10^ (a > 100 A), calculated using Mie 

theory with the 1/3-2/3 approximation (Draine & Malhotra 1993) and dielectric functions from 

Draine & Lee (1984). 

For silicates, we use Mie theory with dielectric constants for astronomical silicate (Draine &: 
Lee 1984).^ 



^We use the "smoothed UV" silicate dielectric function discussed by Weingartner & Draine 
2001a, with modifications to the far-infrared emissivity proposed by Li & Draine 2001a, available at 
http://www.astro.princeton.edu/~draine . 
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Fig. 6. — Optical properties of neutral and ionized PAHs, using C40H16 as an example (see Li & Draine 2001a for 
further information). Also plotted for comparison is the absorption cross section of graphite. Note the difference 
between neutral and ionized PAHs, in particular in the 3-12/im region. In the visible to far UV, the principal 
difference is the location of the visual absorption edge (which is also size-dependent) for ionized PAHs. Insert shows 
the UV absorption for graphite and PAHs. 



4. Excitation of Very Small Grains: Exact Statistical Treatment 

4.1. Assumptions 

We characterize the state of the grain by its vibrational energy E. There are too many 
energy levels to consider individually, so we group them into M + 1 "bins" j = 0, ...,M, where 
the j-th bin is [i^j^min) -E'j,max)) with representative energy Ej = (£'j,min + -E'j,max)/2, and width 
AEj = Ej^yaa.x — Ej^ai\a- The "degeneracy" gj is the number of distinct quantum states included in 
bin J. The procedure used for specifying the bins is described in Appendix B. 

Let Pj be the probability of finding the grain in bin j. The probability vector Pj evolves 
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according to 

lp, = J2TijPj-Y,TjiPi ^ = 0,1,...,M, (14) 

where the transition matrix element T is the probability per unit time for a grain in bin i to make 
a transition to one of the levels in bin /. All of the physics is embodied in the transition matrix T. 

A truly exact treatment of the grain excitation would require knowledge of all the state-to- 
state transition probabilities between energy levels. Given the large number of energy levels to be 
considered, this is utterly infeasible at this epoch. In order to estimate the required state-to-state 
rates, we make only the following three assumptions: 

• The grain absorption cross section Ca_\^s{hi') depends on the photon energy hv, but does not 
depend on the vibrational energy of the grain prior to the absorption. 

• The vibrational modes of the grain can be approximated by harmonic oscillators. 

• The energy of absorbed photons is distributed ergodically among the 3Na — 6 vibrational 
degrees of freedom before any infrared emission takes place.^ 

We now show that with only these three assumptions, together with knowledge of 

• the spectrum of fundamental vibrational modes, 

• the grain absorption section Ca,hs{hi'), and 

• the starlight energy density ue, 

wc can fully determine the emission from interstellar dust grains. Note that the notion of "grain 
temperature" does not appear in what we will refer to as the "exact-statistical treatment" . 



4.2. Upward Transitions 

In Appendix C we show that the rate for upward transitions I ^ it is 



cAE /" 

T^i = |r / Gui{E)C^^,{E)uEdE ioY u < M (15) 

/ — C^UE)uEdE + / C,^,{E)uEdE (16) 



Em -El 



^Since the grain will in general have angular momentum, and the principal values of the moment of inertia tensor 
will generally be nondegenerate, in principle some of the heat can appear in rotational excitation (Lazarian & Draine 
1999a,b). However, since there are only three rotational degrees of freedom, this is at most a small effect, and will 
be neglected. 
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^"'^^^ ^ ^E^AEi ^^'Wi<E<W2 (17) 

- ^^^#1^ fo.^.<.<^3 (18) 

" Xl:^ fa < £ < VF. (19) 

= ioT E <WiOT E>W4^ (20) 

W^l = Eu^min — Ei-^g,^ (21) 

= mill [{Eu ,min -E'i, mill) ) (-^M, max El ,max )] (22) 

W3 = max [(i?,, luin — E/^min), (-E'w.max — -E'/,max)] (23) 

^4 = -E^u.max " -E^Z,min (24) 

Wc = EM,min - -E^Z.min (25) 

and UEdE is the energy density due to photons with energies in [E,E + dE]. Note that eq. (16) 

takes energy absorbed in transitions to levels beyond the highest bin and allocates it to the highest 

bin (M). In practice, we ensure that the highest bin is set high enough that the rate of such 
transitions is negligible, and the population Pm is negligibly small. 

The factor Gui{E) represents the correction for finite bin width; if max(A£';, AEy) <C {Ey_—Ei), 
then eq. (15) becomes 

Tul ^ cC^UE)^^^ iovu<M (26) 



Tmz ^ cC^UE)^^^^ + -ET^-ET r C,UE)uEdE . (27) 
-C'M — -C'Z -C'M — -C/Z J Em -El 



Most previous studies (e.g., Guhathakurta & Draine 1989) have used this approximate form, but 
eq. (15-24) make proper allowance for finite bin width and are to be preferred. Correction for 
finite bin width is important when the treatment is applied to grains with radii a > 50 A, since for 
tractable numbers of bins M < 10^, many bins will necessarily be broad, with A£^„ > kTu- 

For the special case of transitions u—1 u we include "intrabin" absorptions: 

J AEuGu,u-i{E)C^i>s{E)uEdE + J h-——-\c^^^(E)uEdE 

(28) 

but the "intrabin" contribution (second integral in eq. 28) is negligible in all cases of interest. 



-U,Vr-l 



Eu — Eu-i 



4.3. Downward Transitions 



The radiative coupling between levels is determined by the absorption cross section Cabs (^2^) • 
The Einstein A coefficient for radiative decay u —>■ I, averaged over the sublevels in bin u, is (see 
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Appendix C) 



Stt gi AEu 



Gui{E)E^Cai,s{E) 



1 + 



8ttE^ 



ue 



dE for Z < n - 1. 



(29) 



The degeneracies gu and gi are the numbers of energy states in bins u and I, respectively (see §2.3): 



gj = N{Ej^^^) - N{Ej^^i^) « 



dN 
'dE 



AEi 



(30) 



The term proportional to in eq. (29) is the contribution of stimulated emission. For transitions 
n ^ n — 1 to the adjacent bin we have Wi = 0: even very long wavelength emission can contribute 
if C'abs > 0. For this case we must also include the power radiated in "intrabin" transitions within 
bin u: 



-u—l,u 



Stt gu-i 



h^c^ gu E, 
Stt 1 

/l3c2 Eu - Eu- 



AE f^'^ 

] ^ / Gu,u-l{E)E Cabs{E) 

u ~ -C/u-l JO 



1 + 



1 Jo 



1 - 



E 
AE. 



E^CghsiE) 



h?(? 
3 



UE 



1 + 



SttE^ 



Ue 



dE + 
dE 



(31) 



Thus wc sec that wc require only C.ii,^{E), the degeneracies gj, and the starlight spectrum ue to 
completely determine the transition matrix Tjj. We stress that "grain temperature"' does not enter 
the exact-statistical treatment (e.g., eq. 15,29). In §6 below we will discuss "thermal" approxima- 
tions to estimate the downward transition rates T^^. We first discuss an estimate T{E) for the 
"temperature" of a state of specified vibrational energy E. 



5. Energy vs. Temperature 

5.1. "Exact" Mode Spectrum 

Consider a grain with vibrational energy Ey. If the grain is large, we know that the emission 
can be estimated from the grain "temperature" . We can use eq. (2) to define a "temperature" 
T{Eu) such that E{T) = E^, where E{T) is the expectation value for the vibrational energy E 
when in contact with a heat bath at temperature T. 

When there are only a few vibrational quanta in the grain, however, "temperature" is not a 
well-defined concept. As the most extreme example, consider a grain where only the first vibrational 
mode is excited, with exactly one quantum: E^ = fkJi. If the grain has many vibrational degrees of 
freedom, then the temperature T may be quite low (since it is defined in terms of the expectation 
value of the vibrational energy summed over all the modes). However, our grain is vibrating in 
the fundamental vibrational mode and we want to estimate the emission from this mode. One 
estimate for the temperature characterizing the excitation of this mode would be the temperature 
for which the occupation number of this mode is equal to 1. We will see in §6 that when this 
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"temperature" is used in the "thermal approximation" , one obtains the correct rate of emission in 
the 1^0 transition. In fact, we will use this value of the temperature for all energies up to the 
energy of the 20th vibrational mode. For large grains and large energies, however, T is appropriate. 
To encompass these two limits, we take the temperature T„ of vibrational level u to be^ 

kin 2 

= f{Eu) iorEu>hiJ20 (32) 

While the choice of ?7xJ20 as the energy separating these two approximations is arbitrary, we remind 
the reader that the very notion of "temperature" is questionable for such low degrees of excitation. 
Our use of eq. (32) is in fact guided by the fact that the "thermal approximation" of §6, when used 
with this prescription for T^, gives emission rates in fair agreement with the "exact-statistical" 
treatment (29). Figures 7 and 8 show our adopted T{E) for selected PAHs and silicates.^ 



5.2. Debye Mode Spectrum 

When the number of modes is large, the summation in equation (2) contains many terms. In 
this case, for PAHs we replace the sums over the C-C modes by the (continuum) Debye model 
discussed above: 

3 ^ 

^PAH(r) « iVH e^p(^ ./fer) _ 1 + (^c - 2) [kQopf2{T/eop) + 2^9,^/2 (r/e.^)] , (33) 

where /2 is defined by eq. (10), and the summation index j = 1 — 3 runs over the C-H out-of-plane 
bending, in-plane bending, and stretching modes (see §2.1). Similarly, for silicate grains we take 

E,ii{T) « {Na - 2) [2^62/2(^/62) + fee3/3(r/e3)] (34) 

where G2 = 500 K and 63 = 1500 K. We use eq. (33) and (34) for E > NaEi. 



6. Thermal Approximations for Spontaneous Emission 

Given the above assumptions about the vibrational mode spectrum (§2), the assumption that 
the absorption cross section depends only on the photon energy (§3), and the assumption of rapid 
internal vibrational redistribution, the "statistical" approach described above (§4) is "exact" . It is 



^For T = hu}/k\n2 a quantized harmonic oscillator has occupation number equal to 1. 

*The discontinuity in T{E) at E = H0J20 is obviously unphysical, but we note that the discontinuity is only 
appreciable for large grains at very low energies, and these make a negligible contribution to the overall emission 
spectrum. 
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E/hc (cm-') 

Fig. 7. — Temperature T{E) for selected PAH grains with vibrational energy E. For each grain size the dot indicates 
the energy of the lowest vibrational mode. For each grain, T{E) is discontinuous at E = fvjj2o [see eq. (32)]. 

often desirable, however, to use an alternative "thermal" approximation which is less computation- 
ally demanding. 

The thermal approach has been frequently used (Greenberg 1976; Purcell 1976; Aannestad & 
Kenyon 1979; Leger & Puget 1984; Draine & Anderson 1985; Puget et al. 1985; Desert et al. 1986; 
Dwek 1986; Guhathakurta & Draine 1989; Aannestad 1989; Manske & Henning 1998). The only 
difference between the thermal approach and the statistical approach concerns the calculation of 
the downward transition rates which, in contrast to the exact-statistical treatment [eq. (29)], 
uses the notion of "grain temperature" . 

We will consider two different approximations which make use of the thermal approximation, 
and compare results for energy distributions and infrared emission spectra. 
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Fig. 8. — As in Fig. 7, but for silicate grains containing iV atoms. 



6.1. "Thermal-Discrete" Approximation 

In the "thermal-discrete" approximation we use the same energy bins used in the exact sta- 
tistical treatment discussed in §4, and allow the same upward and downward discrete transitions. 
However, instead of using eq.(29) to evaluate the downward transition probabilities T;„, we replace 

91 , 1 

gu AE^exp{hi^/kTu)-l ^ ' 

where T„ is the "temperature" of the upper level. With this replacement, the downward transition 
rate becomes, for / < — 1:^ 



rptd 

^lu — 



87r A£:/ 
^3c2 - El Jw, 



exp(E/A;T„) - 1 



1 + 



8ttE^ 



ue 



dE for < I < It - 1 . (36) 



^Siebenmorgen et al. (1992) discussed both discrete heating and cooling. However, their downward transition 
elements were evaluated by integrating over the initial, rather than final, enthalpy interval, and hence are in error 
by a factor ~ AEu/ AEi, which is large when AEu ^ AEi, as is often the case. Manske & Henning (1998) also 
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For I = u — 1 we include the energy-loss in "intrabin" transitions, as in eq. (31) 

Stt AEi f^^ „ E^C^^,{E) L h^c^ 



-.td 
- u—\,u 



Stt 



j Gu,u-i{E) 
- Eu-1 Jo 



1 + 



ex.p{E/kTu) - 1 
E \ E^C,^,{E) 



Ue 



/^Eu) exp(E/A;T„) - 1 [ ^ttE^ 



1 + 



dE + 

,3 



7UE 



dE . (37) 



We refer to eq. (36,37) as the "thermal-discrete" approximation. The advantage of the "thermal- 
discrete" approximation is that one does not require explicit knowledge of the degeneracies gu and 
gi appearing in equation (29,31). We obtain the temperature from eq. (32).^ 

As discussed in §4.1, the vibrational ground state is in a "bin" with energy Eq = and width 
AEq = 0, for which eq. (36) would give Tq^ = 0. There is really no "correct" way to approximate 
the energy levels as continuous when dealing with the ground state. We simply replace A£^o in eq. 
(36) by AEiJ Thus we do not use eq. (36) for transitions to the ground state, and instead take^ 



ptd 



Stt AEi 1 Z"^"' 
h^AE^K 



E^ 



expiE/kT^) - 1 



1 + 



Ue 



C^bs{E)dE 



(38) 



6.2. "Thermal-Continuous" Cooling Approximation 

In all of the above discussion, the grain is assumed to make discrete transitions to energy levels 
Z < n by emission of single photons - so-called "discrete cooling" . With discrete cooling included, 



discussed discrete cooling. In the present notation their equation (4) becomes 

Comparison with cq. (36) reveals that Manske & Henning's integrand is off by a factor ~ {Eu — 
Ei)/[AEuAEiGui{E)] « {Eu — Ei) / mm.{AEi, AEu) ■ This is a large error for the common case where the pho- 
ton energy is large compared to the energy bin width, E min(A£(, AE^)- 

®Note that the thermal approximation (36) reduces to the statistical treatment (29) if the temperature is taken 
to be 

rp ^ hv/k 

In [1 + {dN/dE)\E„/{dN/dE)\E^-h.] 

but this definition of depends on both u and the emitted photon energy hif - only in the limit of infinite degrees 
of freedom does Tu no longer depend on the value of hv. 

^If this seems rather arbitrary, we remind the reader that the "correct" approach to the problem is the statistical 
treatment given in §4. The "thermal approximation" considered here is inappropriate when discussing low degrees of 
excitation. 

*It can be verified that if we use AEi for the width of bins and 1, and go = gi — 1, then Ti = £i/fcln2 from 
eq. (32) gives the same transition rate Toi as the exact-statistical treatment. 
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the transition matrix Tjj is generally nonzero through the upper triangular portion i > f. 

There are substantial computational advantages (see §8 below) if the cooling of the grains is 
approximated as continuous rather than discrete, so that the only downward transition from a level 
u is to the adjacent level u—1. We refer to this as the "thermal continuous" cooling approximation. 
We take 

Tf^ = if / < -u - 1 (39) 
Using (36), eq.(El), and eq. (E2), the transition rate (40) from levels u > 2 can be approximated 



Because of the computational advantages (see §8.2), the continuous cooling approximation has 
been often used in previous studies (e.g., Guhathakurta &; Draine 1989). While clearly appropriate 
for large systems with many degrees of freedom (for which kT <^ E and emission of one photon does 
not result in a large reduction in temperature), validity of the continuous cooling approximation 
for grains with < 10^ atoms is certainly not obvious a priori. Manske & Henning (1998) in fact 
argued that the continuous cooling approximation could lead to significant error. We will examine 
the accuracy of the continuous-cooling approximation in §10 below, where we will show that it is 
in fact surprisingly accurate even for grains with as few as ~ 30 atoms. 



7. Cooling Time 



The absorption cross section Cabs {E) determines both the rate of heating by starlight and the 
rate of cooling by infrared emission. It is useful to define a "radiative cooling time" 



E„ 



Y^Ku T/w(£^« - El) 
Stt 1 1 rEu 



h^c^ Eu (dN/dE) 



'dN\ 
dE 



dx 



Eu-x 



(42) 
(43) 



where eq. (E2) has been used to obtain eq. (43). In the thermal-discrete approximation, we see from 
eq. (41) that 

1 87r E^C,y,,{E) 



E^ h^c^ 



f 

Jo 



-dE 



for u > 1. 



(44) 



eME/kTu) - 1 

Cooling times for PAH and silicate grains are plotted in Figs. 9 and 10. The interstellar starlight 
radiation field (ISRF) of Mathis, Mezger &; Panagia (1983, hereafter MMP) has been assumed for 
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Fig. 9. — The radiative cooling time r^ad as a function of vibrational energy E for PAH grains, calculated with the 
exact-statistical and thermal-discrete cooling models. The lowest energy E shown for each grain corresponds to the 
lowest vibrationally-excited state hui . The discontinuity in r^ad for the thermal-discrete model is due to the different 
definition of the grain "vibrational temperature" for the lowest 20 excited levels (see §5 and Fig. 7). Note that r^ad 
for the thermal- continuous model is identical to that for the thermal-discrete model [eq. (15) and eq. (36)]. 

ue', for this very dilute radiation field the correction for stimulated emission is entirely negligible 
{h^c^UE/87rE^ < 1 for E/hc > lOcm-i). 

We show Tj-ad evaluated using the "exact-statistical" method, as well as using the "thermal- 
discrete" approximation (rrad foi" the continuous cooling approximation is, by construction, identical 
to that for the "thermal-discrete" approximation). 

As illustrated in Figs. 9-10, TradC-E') increases with grain size; for a given grain size, r^ad de- 
creases rapidly with increasing vibrational energy E. The importance of single photon heating 
is apparent: very small grains emit most of their energy very rapidly (within seconds) following 
absorption of an energetic photon and spend most of their time in or close to the ground vibrational 
state. 
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Fig. 10. — Same as Fig. 9 but for silicate grains. 

The substantial difference between the results of the statistical model and the thermal models 
at low E is mainly due to the treatment of transitions to the ground state, which are problematic 
in the thermal approximation (see §6.1). The jagged behavior of Tj-ad in the exact-statistical model 
result from the discreteness of the lowest energy states (see eq. 4) which affects the degeneracies 
(see eq. 29). The discontinuity in r^ad computed in the thermal approximation (Figures 9,10) is due 
to the discontinuity in the prescription (32) for the grain "vibrational temperature" at E = huJ2o 
(see Figs. 7,8). Note that with the present prescription for the temperature, our estimates for 
Trad in the thermal approximation are in reasonable agreement with r^ad calculated in the "exact- 
statistical" method at the lowest energies. For C4000H1000 in Figure 9, and a = 20 A in Figure 10, 
the discontinuity in Trad at ~ 3 x 10^ cm~^ is due to a transition from direct computation of the 
density of states to use of a continuous Debye model. 



The mean time Tabs 

between photon absorptions for a grain is given by 



(45) 
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In Figure 11 we show r^bs for grains in the MMP ISRF. The mean absorbed photon energy (/iJ^)abs 
is 

{hy)shs = Tabs / Cshs{hv)cui,dv . (46) 
Jo 

Ionized PAHs with A^'c < 10^ in the MMP radiation field have {hu)shs ~ 5.2 eV; neutral PAHs 
have a somewhat higher value of {hi')abs due to their different absorption cross section (see Figure 
6). 



IC 



106 



105 



10" 



J3 
(S 



10= 



,^102 

10 

1 

0.1 



a (A) 

30 50 



10 




Fig. 11. — Radiative cooling time Trad for a grain containing E = 5.2 eV of vibrational energy, and mean time Tabs 
between photon absorptions for the MMP ISRF, as a function of grain size, for carbonaceous and silicate grains. The 
decrease in Trad as A^'c increases above ~ 5 x 10'' is due to the assumed transition from PAH-like to graphitic optical 
properties at a « 50 A. 



Also plotted in Figure 11 is the cooling time Trad for grains containing 5.2 eV of thermal energy. 
The critical size Agph (and radius Ogph) for single-photon heating is determined by the condition 
Trad(A^sph) = Tabs(Asph)- For carbonaccous grains, we find 



A^Csph « 5 X 10^ 



-0.60 



, asph-lOOA 



^iMMP 



-0.20 



(47) 
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For grains with N < ATgpjj atoms, single-photon heating effects are very important: the radiative 
cooHng time is short compared to the interval between photon absorptions, and the grains spend 
most of the time between photon absorptions at energies E <^ (/ii^)abs- Grains with N S> Nsph heat 
up until their energy content is large compared to the energies of starlight photons, and individual 
photon absorptions produce small fractional changes in the thermal energy content. For such large 
grains, one may therefore estimate the average energy content E by requiring a steady balance 
between cooling and heating: 

E f 

= j dv cUyCabsihv) . (48) 



8. Solution Method 
8.1. Discrete Cooling 

If we define the diagonal elements of T to be 

T,, = -5^T,, (49) 

then equation (14) becomes 

M 

^ TijPj = for z = 0, M . (50) 

Combining this with the normalization condition X^^q — obtain a set of M linear equations 
for the first M elements of Pj: 

M-l 

^ (Ti, -TiM)Pj = -TiM fori = 0,...,M-l . (51) 

i=o 

The remaining undetermined element Pm is obtained by 

M-l 

Pm = - (Tmm)"^ J2 ^^j^j ■ (^2) 

j=0 

For small M we can directly solve eq. (51) for Pj using Gaussian elimination. When M > 10^, 
however, direct solution, requiring O(M^) operations, is both slow and numerically ill-behaved. 
We therefore resort to iterative techniques, and have tried both the bi-conjugate gradient (BiCG) 
method (Press et al. 1992), and the stabilized bi-conjugate gradient (BiCGstab) method (see 
Sleijpen & Fokkema 1993 and references therein). We found the BiCG method to be quite efficient 
even without preconditioning, whereas the convergence rate of the BiCGstab method is dependent 
on the choice of preconditioner. 
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Ideally, the adopted highest energy bin Em for a given grain would be as large as possible, 

but large values of M entail heavy memory and cpu requirements. In our calculations, we first set 
Em = 13.6 eV and solve eq. (51) for Pj. If Pm < 10~^^ then we use this Em and solution Pj; 
otherwise we increase Em by a factor 1.5 and solve again, repeating until we obtain a solution with 
Pm < 10-^1 

SmaU (a < 20A) grains exposed to the interstellar starlight background are able to cool 
completely betweeen photon absorptions, and Em = 13.6 eV is sufficient. For larger grains, the 
interval between photon absorptions is shorter and the cooUng time at given E is longer (see 
Figs. 9 and 10), so that the grain generally does not cool off to the ground state between photon 
absorptions, and Em must be higher than 13.6eV. For a grain with a > lOOA only modest cooling 
takes place between photon absorptions and Pj becomes strongly peaked around a steady state 
"equilibrium temperature" . 



8.2. Continuous Cooling 

For the continous cooling approximation it is straightforward to solve for the steady-state 
solution vector Pj. As shown by Guhathakurta & Draine (1989), if we define Xj = Pj/Pq, then 
Xq = 1 and the remaining Xj may be obtained recursively: 

^ i-i M 

i=0 u=j 

The recursion relation Bj^i = Tj^i+Bj+i^i can be used to generate the Bj^i in O(M^) operations, and 
the Xj can then be computed in O(M^) operations. Once the Xj are determined for j = 0, M, 
one need only renormalize: 

Since M > 100 is required to have suitably narrow energy bins (we often use M « 500), this 0(M) 
procedure is much faster than the BiCG iterative method required if discrete cooling is included. 



9. Infrared Emission Spectrum 

Having obtained the steady-state energy distribution Pj, we wish to calculate the resulting 
infrared emission spectrum. Let Fi,di' be the time-averaged power per steradian per grain radiated 
into frequency interval [u, v -|- dv\ . Then the "exact statistical" treatment gives 



2/ji/4 



u-1 



J^Pu X -^EuGuiihu) + J^'Pu ( 1 - 



1=0 



9u 



hv 
AE. 



1 + -^UE 



(55) 
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The single-primed summation - the contribution of "interbin" transitions to the emission - is limited 
to levels u with Eu^max > hu-^ The double-primed summation - the contribution of "intrabin" 
transitions - is limited to levels with AE^ > hv. 



In the thermal approximation, equation (55) becomes (see Appendix D) 



eicp{hi^/kTu) - 1 

where the single-primed summation is now limited to levels u with Eu > hu. 



(56) 



10. Results 

Wc have calculated the (vibrational) energy probability distribution Pj for PAHs excited by 
the MMP ISRF using the above-described methods: the exact- statistical model (§4); the thermal- 
discrete model (§6.1); and the thermal- continuous model (§6.2). 

In Figure 12 we present the cumulative energy probability distributions for selected PAHs 
obtained from the exact-statistical model, the thermal-discrete model, and the thermal-continuous 
model. Figure 13 shows the differential probability distributions as a function of grain energy E 
for these same grain sizes. 

The inset in Figure 13 shows the probability of being in the vibrational ground state as a 
function of grain size, calculated for the exact-statistical model and the thermal-discrete model. As 

already illustrated in Figure 12 and Figure 13, the probability of being in the ground state is very 
large for small grains: for example, for the MMP radiation field, grains with N < 4000 spend most 
of their time at £^ = [the exact-statistical model gives Pq = 0.975 for A'^ = 400, and Pq = 0.316 for 
N = 4000]. The sharp drop at 13.6 eV {E/hc = 1.1 x 10^ cm-^) is due to the radiation field cutoff 
at 912 A and to the fact that multiphoton events are rare. The jagged structure at low energies in 
Figure 13 is due to the quantization of the lowest energy states. These structures arc less prominent 
in the thermal models because the statistical model involves the degeneracies of each "bin" (see 
eq. 30) which depend on the detailed distribution of energy states (see, e.g., Figure 4). 

The energy distributions P{E) found using the thermal-discrete model and the thermal- 
continuous model are both in good overall agreement with the results of the exact-statistical cal- 
culation. 

It is not surprising that all three models are in good agreement at high energies since the 
transition rate from a very high energy state to the ground state is negligibly small, and the 



®Note that only a few terms from Y^'^ have G„; ^ - energy levels k <l <m where i?u,min — hv£ [-Bfc.min, £fe,max] 
and Eu,max — hv € [Em,min, Em,max]- Therefore eq. (55) requires only 0{M) [rather than 0(M^)] operations per 
frequency v. 
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Fig. 12. — The cumulative energy probability distributions for PAHs computed using the exact-statistical model, 
the thermal-discrete model, and the thermal-continuous model. Note that the lowest energy state {E = 0), not shown 
here, has P(E > 0) = 1. 

number of excited degrees of freedom is large enough that a thermal treatment is appropriate. 
The agreement becomes better with increasing grain size because the lowest vibrational states (i.e., 
the first few excited states) become closer and approach a "continuum" , and the gap between the 
ground state and the first excited state becomes smaller. 

The resulting IR emission spectra are displayed in Figs. 14 and 15. As expected (from the 
energy distribution functions plotted in Figs. 12 and 13), the IR spectra of the thermal-discrete 
model are almost identical to those of the exact-statistical model, even for grains as small as 
Nc = 24. The thermal-continuous cooling model also results in spectra which are very close to 
those computed using the exact-statistical model. The discrepancies are mainly at long wavelengths; 
this can be understood from the energy probability functions and would not affect the interstellar 
IR emission spectrum modeling. 

The statistical model spectra are generally very close to those of the thermal models (Figs. 14 
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Fig. 13. — The vibrational energy distribution for excited vibrational states of ionized PAHs in the MMP ISRF, 
computed using the exact-statistical model, the thermal-discrete model, and the thermal- continuous model. The inset 
shows the probability Pq of being in the ground state. 

and 15), except for the sawtooth features at long wavelengths. These sawtooth features are due to 
our treatment of transitions from the lower excited energy bins to the ground state and first few 
excited states. However, we stress that the overall spectra are quite similar for all three models - 
the differences involve only a negligible fraction of the total emission for any given grain size, and 
thus would be unimportant when modeling the overall interstellar IR emission spectrum. 

11. Discussion 
11.1. Band Ratios 

From Figs. 14 and 15 it is evident that the relative strength of the different PAH emission 
bands (3.3, 6.2, 7.7, 8.6, 11.3/xm) is a strong function of the PAH size: small PAHs radiate strongly 
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Fig. 14. — IR emissivities (per C atom) for selected ionized PAHs in the MMP ISRF calculated using the exact- 
statistical and thermal-discrete models. 

at 6.2 and 7.7/im, while larger PAHs emit most of their power at increasingly long wavelengths. 

The 6.2, 7.7, and ll.S^m features are prominent in many astronomical spectra. Let -P(Ao) 
be the power radiated in the feature with central wavelength Aq: for a Drude profile, -P(Ao) = 
(7r/2)(AP\)o(FWHM)AQ, where (APa)o is the peak contribution of the feature to the emission, 
and (FWHM)ao is the full-width at half-maximum. In Figure 16 we show how the relative 
strengths of these emission bands vary depending upon the size and charge state of the PAHs, 
and upon the starlight intensity ue-, characterized by x oi" G^}^ Neutral PAHs have large ra- 
tios of P(11.3^m)/P(7.7/im), while PAH ions have much smaller values. For both neutrals and 
ions, there is a regular progression of P(6.2/im)/P(7.7/um) to smaller values with increasing A'^c- 
For a given PAH neutral or ion, the emission band ratios are essentially independent of x for 



^'-'x is the intensity at 1000 A relative to the estimate of Habing (1968) ; Go is the ratio of the 6-13.6 eV energy density 
to the estimate of Habing (1968). The MMP spectrum has Go = 0.923x; a 3 x 10* K blackbody has Go = 0.602x. 
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Fig. 15. — IR emissivities (per C atom) for ionized PAHs of calculated using the thermal-discrete and thermal- 
continuous cooling models. 

[see eq. (47)]; for larger values of x, multiphoton heating effects begin to 
be evident. For x = 10^ i multiphoton effects are evident even for Nq = 10^ [consistent with eq. 
(47)]. 

In Figure 16 we plot the observed ratios of P(11.3/xm)/P(7.7;um) and P(6.2/im)/P(7.7;um) for 
(1) the diffuse ISM (Onaka et al. 1996); (2) the average spectrum for a sample of 28 spiral galaxies 
(Helou et al. 2000); (3) the reflection nebula at the surface of the p Oph molecular cloud (Boulanger 
et al. 1996); (4) the reflection nebula vdB 133, illuminated by F5Iab and B7II stars (Uchida et 
al. 1998); (5) the reflection nebula/photodissociation region (PDR) NGC 7023, illuminated by a 
B2.5V star (Cesarsky et al. 1996a); (6) the M17 PDR (Cesarsky et al. 1996b); (7) the starburst 
galaxies M82 and NGC 253 (Stiirm et al. 2000); (8) the Seyfert 2 galaxy in Circinus (Stiirm et al. 
2000); and (9) the quiescent molecular cloud SMC Bl^l in the Small Magellanic Cloud (Reach et 
al. 2000). With the exception of SMC Bl#l (which we discuss below), the observed band ratios 
fall between the neutral and ionized PAH tracks in Figure 16, suggesting that the observed spectra 
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Fig. 16. — Relative strengths of three PAH emission features (6.2 and 7.7/im C-C stretching modes, and 11.3/im 
C-H out-of-plane bending mode) for neutral and ionized PAHs, labelled by the number of carbon atoms. Results are 
shown for grains illuminated by the MMP spectrum with x ~ 1-23 and 123, and for a 3 x 10^ K blackbody spectrum 
(cutoff at 912 A) with x = 10^ 10*, 10^, and lO'^. Also shown (stars) are observed band ratios for the diffuse ISM, 
the reflection nebulae NGC 7023 and vdB 133, the M17 PDR, the p Oph molecular cloud; the starburst galaxies 
M82 and NGC 253 (Stiirm et al. 2000); the Seyfert 2 galaxy in Circinus; the average of 28 normal galaxies; and a 
quiescent molecular cloud in the SMC. (See text for citations.) The dotted line connecting the neutral and ionized 
points for A'^c = 16 (in a PDR) illustrates the variation in the relative strengths of the 6.2, 7.7, and 11.3/im features 
in a PDR as the PAH ion fraction varies from to 1. 

can be reproduced by an appropriate mixture of PAH sizes, with an appropriate ionized fraction. 
The observed band ratios for the Milky Way ISM, p Oph, 28 normal galaxies, and the reflection 
nebulae vdB 133 and NGC 7023 could be reproduced with a relatively low ionized fraction, but 
larger ionized fractions are required for the M17 PDR, M82, NGC 253, and the Circinus galaxy. 

11.2. SMC Bl#l 

The quiescent molecular cloud SMC Bl^^l observed by Reach et al. (2000) falls outside the 
region bounded by the neutral and ionized tracks in Figure 16. Reach et al. point out that SMC 
Bl^l has a 11.3/7.7 band ratio which differs significantly from that for Milky Way objects, includ- 
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• CNM (MMP ISRF. n„=30cm-=. T=10bK, nyn„= 1 .5x 10"') 
O same as CNM but with lOOxMMP ISRF. n„=3000cm-= 

• WIM (MMP ISRF, n„=0.1cm-=, T=8000K, n,/n„=0.99) 
O same as VflM but with lOOxMMP ISRF, n„=10cm-" 

• PDR (T.=3xl0-K, x=10^ n„=10=cm-^ T=1000K, n,/n„=3x lO-) 
□ same as PDR but with x=lO'' n„=10»cm-= 



SMC 




0.06 
0.1 



0.5 



0.2 0.3 
P(6.2/xm)/P(7.7yu,m) 

Fig. 17. — Relative strengths of three PAH emission features (6.2 and 7.7/im C-C stretching modes, and 11.3/im 
C-H out-of-plane bending mode); PAHs are labelled by the number of carbon atoms. Results are shown for 3 charging 
conditions: CNM, WIM, and PDR (see text). Circles: grains in CNM charging conditions, for Go ~ 1.14 and 114; 
diamonds: grains in WIM charging conditions, for Go ~ 1.14 and 114; squares: grains in PDR conditions, for x = 10^, 
10*, 10^, and lO''. Stars represent observations as in Fig. 16. 



ing the p Oph reflection nebula (which one might expect SMC Bl#l to resemble). Reach et al. 
suggest that the enhanced emission in the 11.3/im C-H bending mode might be due to increased 
hydrogenation of the PAHs in the SMC, where the gas phase H/C ratio is a factor of 10 higher 
than in the Milky Way. Given that Milky Way PAHs with >25 C atoms are expected to be essen- 
tially fully hydrogenated (Tielens et al. 1987; Allamandola, Tielens, &: Barker 1989; Allain, Leach, 
& Sedlmayr 1996), it is not clear that greater hydrogenation would be expected in the SMC. As 
discussed in Paper II, our model calculations assume that interstellar neutral and ionized PAHs 
have the 11.3/7.7;um band strength ratio reduced by a factor of 2 relative to PAHs which have been 
studied in the lab - enhancement of the 7.7pm feature appears to be necessary to reproduce the 
low 11.3/7.7/im band ratios observed in objects such as NGC 253 or M82, and (with the possible 
exception of SMC Bl^^l) the reduced 11.3/7.7;um band strength ratio appears to be consistent 
with astronomical observations. 

While SMC Bl^^l appears to fall "outside" the "allowed" region in Figures 16-18, we show 
elsewhere (Li & Draine 2001b) that the present model can reproduce the observed line ratios by 
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0.1 0.2 0.3 0.4 

P(6.2//m)/P(7.7yLim) 

Fig. 18. — Same as Fig. 16 but for a weighted sum of emission from grains in the cold neutral medium (CNM), 
grains in the warm neutral medium (WNM), and grains in the warm ionized medium (WIM), heated by the MMP 
estimate of the ISRF. Band ratios within the shaded region can be obtained by appropriate mixtures of PAHs with 
20 < A^c < 1000. Observed band ratios for the diffuse ISM (Onaka et al. 1996; Mattila et al. 1996) and the average 
spectrum of 28 normal spiral galaxies (Helou et al. 2000) fall within this allowed region. 

adding emission from a population of larger PAHs - with Nq ~ 3000 - which radiate strongly in 
the ll-13//m region, but produce little emission in the 6-8/im region. These larger PAHs would 
radiate strongly at 15-50/im; SIRTF will be able to test this hypothesis. 

11.3. PAH Ionization 

The fraction of PAHs which are charged will depend upon the electron density Ug, temperature 
T, and the starlight intensity. We consider three different charging environments: 

• "Cold neutral medium" (CNM; = 0.045 cm'^ T = 100 K, and the MMP ISRF spectrum 
with X = 1-23, Go = 1.14) and with both radiation field and gas density increased by 10^; 
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• "Warm ionized medium" (WIM; rig = 0.1 cm'^, T = 8000 K, and the MMP ISRF with 
X = 1-23), and with both radiation field and gas density increased by 10^; 

• Photodissociation region conditions (FDR; Ue/riB. = 3 x 10~^, T = 10^ K, and a 3 x lO^K 
blackbody cut off at 912 A with x/nu = 0.1 cm^) for x = 10^, 10^, 10^, and 10^. While PDRs 

come in a range of conditions, recent observations and modeling suggest that at least some 
PDRs include an extended region with T ~ 10'^ K, and with x/'^H ~ 0.1 cm^, as in a recent 
model for the NGC 2023 PDR (Draine k Bertoldi 2000). 

For each set of conditions, we consider PAHs with Nc > 16 C atoms, and include emission from 
both neutral and ionized forms, with the neutral fraction calculated assuming a balance between 

collisional charging by electrons and photodetachment and photoionization (Wcingartncr & Draine 
2001b). The charging of small grains is largely a function of the quantity x^T/ue (Bakes &: Tielens 
1994; Weingartner & Draine 2001b). For the CNM, WIM, and PDR conditions considered here, 
xVr/ne = 270, 1100, and 10^ cm^K^/^, so we expect grains of a given size to be increasingly 
positively charged as we move from CNM to WIM to PDR. 

The "diffuse ISM" and "28 Galaxies" points fall between the CNM and WIM tracks in Figure 
16, implying that one could reproduce the overall emission from normal spiral galaxies by a weighted 
sum of interstellar conditions ranging from CNM to WIM, with a mixture of grain sizes such that 
the 6.2 - 11.3/Ltm emission is dominated by grains with 20 < Nq < 500. 

The observed flux ratios for NGC 7023 and M17 fall between the "WIM" and "PDR" tracks 

in Figure 17. The observed flux ratios for the starburst galaxies M82 and NGC 253, as well as the 
Circinus Scyfcrt 2 galaxy, fall close to the "PDR" tracks, consistent with the expectation that the 
spectra of these galaxies includes a substantial contribution from PDRs. 

It is important to note that the "PDR" models in Figure 17 have charging conditions such 
that a substantial fraction of the small grains are still neutral: 22% of the Nq = 100 grains are 
neutral. Larger charged fractions - and smaller values of P(11.3/Lim)/P(7.7/Ltm) - can be achieved 
by increasing the value of x/i^H- For example, if x/'^-H is increased by a factor 2, only 10% of the 
Nc = 100 grains are neutral. 

In Figure 18 we show the emission expected from grains of various sizes, with 43% in the 
CNM, 14% in the WIM, and 43% in the WNM {ue = 0.03 cm'^ T = 6000 K, MMP ISRF, 
xVf/ue = 3200). The shaded region shows the flux ratios which could be obtained by grain 
mixtures containing only grains 20 < Nq < 1350 (from Figure 15 it is apparent that grains with 
Nq > 10^ are not likely to contribute significant emission at A < 11.3/im). The observed flux ratios 
for the diffuse ISM in the Milky Way, and for the average of 28 spiral galaxies, both fall in the 
shaded region, suggesting that with an appropriate size distribution, one might be able to reproduce 
the observed emission. The question of size distributions will be addressed in a future paper (Li & 
Draine 2001a), where it will be shown that the observed infrared emission can be approximately 
reproduced by grain size distributions similar to those invoked to explain the observed microwave 



35 



emission from interstellar dust (Draine &; Lazarian 1998a,b). 

12. Summary 

We have presented a method for calculating the infrared emission from dust grains, including 
very small PAH molecules or silicate clusters. The principal results of this paper are the following: 

1. The vibrational mode spectrum and vibrational density of states of very small grains are 
discussed in §2, where we explicitly consider both polycyclic aromatic hydrocarbons and 
silicate grains. 

2. We use the adopted vibrational mode spectrum to calculate the specific heat of carbonaceous 
and silicate grains; the results are in agreement with experimental specific heats for graphite 
and amorphous silicates (Figure 2). 

3. Using estimates for the absorption cross sections of carbonaceous and silicate grains, we obtain 
(in §4) radiative transition rates for emission from vibrationally-excited states of carbonaceous 
and silicate grains. We refer to this as the "exact-statistical" treatment. 

4. We discuss (in §6) how the "thermal approximation" can be used to estimate rates for sponta- 
neous emission. With a judicious estimate for the "temperature" , the thermal approximation 
provides transition rates which are close to those obtained from the "exact-statistical" treat- 
ment. Therefore level populations estimated using the thermal approximation are close to 
the actual level populations given by the "exact-statistical" treatment. 

5. We discuss the "continuous cooling" approximation for modeling the deexcitation of the grain. 
We show that level populations and emission spectra computed using this approximation are 
sufficiently accurate for most astrophysical applications, even for grains containing as few as 
~ 30 atoms. 

6. Methods for numerical solution for the vibrational energy distribution arc presented in §8, 
and the numerical advantages of the continuous cooling approximation are stressed. 

7. Relative strengths of the 6.2, 7.7, and 11.3/xm features depend on the grain size, on the 
charging conditions, and on the starlight intensity. Observed 6.2/7.7 and 11.3/7.7 band 
ratios for the diffuse ISM in the Milky Way, for the p Oph and vdB133 reflection nebulae, 
and for the 28 "normal" galaxy sample of Helou ct al. (2000), appear consistent with grains 
in the diffuse ISM heated by the average starlight background. The band ratios observed 
for the starburst galaxies M82 and NGC 253, for the Seyfert 2 "Circinus galaxy", for the 
M17 PDR, and for the NGC 7023 FDR appear to require conditions close to our "FDR" 
conditions, with x/nH~0.05 — 0.1 cm^. 



36 



8. The unusual band strengths observed for the quiescent molecular cloud SMC Bl#l in the 
SMC (Reach et al. 2000) cannot be reproduced by a mixture of small Nq < 10^ PAHs. Either 
the PAHs in this cloud differ from Milky Way PAHs in their properties, or have a different 
size distribution. 

The present paper has concentrated on the statistical mechanics of the stochastic heating process, 
and on the resulting emission spectra as a function of grain size. In a separate paper (Li & Draine 
2001a) we calculate the emission from grain models with realistic size distributions, and show that 
the observed infrared emission from interstellar dust can be reproduced by a grain model with 
~15% of the interstellar carbon in PAH particles containing less than 10^ C atoms. 

We thank L.J. Allamandola, J.M. Greenberg, and J.C. Weingartner for helpful discussions. We 
thank D.M. Hudgins and L.J. Allamandola for providing us with their PAH database. We thank 
R.H. Lupton for availability of the SM plotting package. This research was supported in part by 
NSF grant AST-9619429, and NASA grant NAG5-7030. 



A. PAH Geometry 

1/2 

Small PAHs are expected to be planar, with D oc NJ , but large PAHs are expected to be 

1/3 . 

approximately spherical, with D oc N^' . The size at which the transition from planar to spherical 
geometry occurs is uncertain. We will assume that PAHs have a planar geometry up to A^c = 54 
(e.g., the pericondensed species circumcoronene C54H18). We assume that beyond A^c = 54 PAHs 

consist of multiple sheets. Beyond Nq = 102 (e.g., C54H18 plus one coronene C24H12 on either side) 

1/3 

we assume spherical geometry, with D oc . 



B. Energy Bins 

We specify the energy bins so that the total number M of bins is manageable (300 < M < 
1500), the number of vibrational states per bin is at least 2 for each of bins 3-10, the bin widths vary 
smoothly for j > 10, and the last bin is at a high enough energy that its population Pm < 10~^^ 
(see §8.1 below). 

Let ujj,j = 1, ...,3(A^a — 2) be the fundamental vibrational modes (see §2), in order of increasing 
frequency. Bin is the ground state: -Eo,min = -E'O.max = Eq = 0, to which we assign a "width" 
AEq = 0. Bins 1 and 2 contain the first two vibrational modes, and bins 3 through 10 each contain 
two normal modes: 

3 1 

El,mm = -jhjJl - 2^^2 (Bl) 
-E^j.max = -Ej+l,min = ^ (^j + ^j+l) for j = 1, 2 (B2) 
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-Ej,max = -Ej+l,min = (^2j-2 + ^2j-l) for 3 < j < 10 (B3) 

For j = 11, ...,K we use uniform intervals AE, and for j > K uniform intervals in ln£^: 

^j,max = £^10,max + (j - 10)(£;io,me« - ^9,max) for 10 < J < if (B4) 

^j,max = £^K,max f Z'^''"" V for K < j < M , (B5) 

\-C/_fs:-l,max/ 

if is chosen so that the total number of bins ~ 500. 



C. Transition Matrix 

Consider energy "bins" u and I, with E^ > Ei. Suppose bin I consists of many subbins 
a = 1, ...,Nci, each of width AE^ = AEi/Na, and bin u consists of many subbinss (3 = 1, ■■■^Nj}, 
each with width AEp = AE^/Np. The upward Einstein A coefficient from subbin a to subbin (3 is 

Apa = -f^^EpC^^,iUap) , (CI) 

and the spontaneous decay rate from subbin f3 to subbin a is 

_ {dN/dE)^ 87r,.lg 

- {dN/dE)p hc^ ^E-C^^^M , (C2) 

where Uafs = {Ep — Ea)/h, and dN/dE is the vibrational density of states. We now assume that 
{dN/dE)a ~ {dN/dE)i, and (dN/dE) ~ {dN/dE)^. 

To obtain the effective bin-to-bin transition rate T„;, we equate the total rate of energy ab- 
sorption to T„/(£'u — El). Thus, averaged over subbins a, the transition rate from Z — >■ n is 

Tui = TT^if EE^abs(^^)^EAi?;3 (C3) 
E^-EiN^^y 



1 rEi ,max rEu,raa.x 

I Ih I U I2j1 Trim " J-^ii mm 



The integrand depends only on E = y — x. One can show that 
where Wi, W4, and Gui{E) are defined in equations (19-24) 



T„/ = -T |r / dEGui{E)C,^,iE)uE , (C5) 
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Similarly, averaged over subbins /?, the energy-weighted transition rate from -u — ^ Z is 



^^^EE^-M^.-^^) (C6) 

{dN/dE)i Svr ^Ep Wf3^ .p7^ 



; / d£;G„K^)i?'Cabs(^) • (C9) 

■u — m Jw^ 



91 Stt 



gu h^c^ Eu - El 

In the case of transitions to adjacent bins, we augment T„^„_i by the rate of energy absorption 
in "intrabin" transitions divided by E^ — E^-i, and T„_i^u by the power radiated in intrabin 
transition, divided by E^ — -E'u-i (see eq. 28, 31). 



D. Emission Spectrum 

Fjy, the power radiated per grain per unit solid angle per Tinit frequency, can be evaluated by 
noting from eq. (29) that the contribution from transitions u ^ I is just 



5F,du = P^^-l^Gui{E)E^C^^,{E) 

gu h-^c^ 



■ A3 ■ 



dE 



(Dl) 



If the width of bin u exceeds hi', then there is an additional contribution from intrabin transitions 
originating in the fraction (1 — hi^/AE) of the subbins which are more than hu cibove Ey^ min- 

Using 

eq. (C2) with {dN/dE)a ~ {dN/dE)p, the intrabin contribution to F^dv is 



SF^du 



Pull- 



= Pu 1 



hu \ E dA 



AEuJ AiidE 
hu \ E 87ri/2 



AE„ / 4tt hc^ 



dE 

Svr 



dE 



(D2) 
(D3) 



where dA/dE is obtained by dividing eq. (C2) by AEa- With E = hv, it is straightforward to use 
eq. (D1-D3) to obtain eq. (55). 



E. Lemmas 



From the definition (19) of Guiihv) one can show that 



u-l 



AEuAEi Guiihu) = min {hu, AE^) for hu < Eu,ram - Ei^, 



1=0 



(El) 
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u-l 

Y,^Ei{E^^Ei)G^l{E)^E , (E2) 

1=0 

which are used to obtain eq. (40) and (43). 
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